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Abstract 



We present an approach for polarimetric Synthetic Aperture Radar 
(SAR) image region boundary detection based on the use of B-Spline 
active contours and a new model for polarimetric SAR data: the Qp 
distribution. In order to detect the boundary of a region, initial B- 
Spline curves are specified, either automatically or manually, and the 
proposed algorithm uses a deformable contours technique to find the 
boundary. In doing this, the parameters of the polarimetric Qp model 
for the data are estimated, in order to find the transition points be- 
tween the region being segmented and the surrounding area. This is 
a local algorithm since it works only on the region to be segmented. 
Results of its performance are presented. 

1 Introduction 

Synthetic Aperture Radar (SAR) data have proven their importance in a 
number of applications, among which can be mentioned forestry, agriculture, 
analysis of geological and geomorphological features, thematic map updating, 
ocean oil spills, marine climatology, ice monitoring and deforestation (see, 
among others, 0, EJ [TO]). 

*A. C. Frery is with the Instituto de Computagao, Universidade Federal de Alagoas, 
Brazil. M. E. Mejail, J. Gambini and J. Jacobo-Berlles are with the Departamento de 
Computation, Universidad de Buenos Aires, Argentina. 



Due to the active nature of the sensor, such images can be obtained at 
any time of the day, since the illumination source is carried by the sensing 
device. SAR sensors work on the microwaves spectrum, so they are almost 
immune to adverse weather conditions and they are able to penetrate, to 
some extent, the surface of certain targets. The first civilian SAR satellite 
was launched in 1978, and it was followed by a constellation of other similar 
sensors, mostly devoted to specific applications and in all cases operated at 
a single frequency and polarization. 

The Shuttle Imaging Radar-C/X-band SAR (SIR-C/XSAR), launched in 
1994, could be operated simultaneously at three frequencies, with two of 
them able to transmit and receive at both horizontal and vertical polariza- 
tion. This polarimetric capability provides a more complete description of 
the target [16]. 

Polarimetric images are multiple complex- valued data sets requiring, thus, 
specialized models and algorithms. Two main venues of research are followed 
for such description: electromagnetic modeling (see, for instance, |37j and the 
references therein), and statistical laws. We follow the second one. 

Many techniques have been proposed for feature extraction and area clas- 
sification in polarimetric SAR imagery. Migliaccio et al. |32j present a study 
on sea oil spill observation by means of polarimetric SAR data, based on the 
use of a constant false alarm rate filter. Some techniques (see, for instance, 
[38J) treat each channel individually and then fuse the results, but this ap- 
proach does not exploit all the information these images convey. Goudail et 
al. [23] present a framework for designing algorithms that can solve detection, 
location and segmentation in polarimetric SAR images. In these works the 
authors propose a definition of the contrast between regions with different 
polarimetric properties. Horta et al. [26J perform polarimetric SAR image 
classification using the EM algorithm. 

Conradsen et al. [I] model polarimetric SAR data with the complex 
Wishart distribution for edge detection. Under that model, Davison et al. [6] 
perform classification by maximum likelihood, while Ferro-Famil et al [2H] 
present segmentation results. Other distributions used for polarimetric SAR 
image segmentation are the K p [B] and the Q% laws [12J, [13] . 

Our approach is also based on the statistical description of the available 
information, which is in the form of a matrix in every pixel. An edge is 
an ideal curve that divides two areas with different polarimetric scattering 
mechanisms, which yield different statistical properties in at least one of the 
available components; in our case, the roughness will be used as discrimina- 
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tory feature. This feature will be described as a real-valued parameter that 
indexes a new distribution for polarimetric SAR data, namely, the Harmonic 
G law (denoted Gp)- 

To the knowledge of the authors of this paper, the literature shows no 
results of combining polarimetric statistical properties with B-Spline based 
deformable contour methods. Our method has the ability of operating on 
regions instead of over the whole image, which is a considerable virtue given 
the complexity and size of polarimetric SAR images. 

The B-Spline approach has been widely used in curve representation for 
boundary detection [37], among other applications. Contours formulated by 
means of B-Splines allow local control of the curve, have local representation, 
require few parameters and are intrinsically smooth. A method for bound- 
ary estimation in noisy images based on B-Spline deformable contours, the 
Minimun Length Criterion and the Gaussian distribution is described in [11] . 
Gambini et al. (TTJ [18] developed techniques for boundary detection in uni- 
variate amplitude SAR imagery using the Q\ distribution. 

The technique proposed in this work is based on B-Spline boundary fit- 
ting, as proposed by Blake and Isard [2j, but tailored to the properties of 
polarimetric SAR imagery by means of the polarimetric Gp distribution as 
a general data model. The polarimetric Gp model was recently developed, 
and presents an attractive choice for polarimetric SAR data segmentation. 

This proposal for boundary extraction begins with the manual or auto- 
matic specification of initial regions of interest, determined by control points 
which generate a B-Spline curve. Then, a series of radial segments are drawn 
on the image, and image data around them are extracted. For each segment, 
the transition point, that is, the point belonging to the region boundary, 
is determined by parameter estimation from the data under the Gp model. 
Then, for each region, the contour sought is given by the B-Spline curve that 
fits these transition points. 

We apply this technique to simulated data in order to quantitatively 
assess its performance, and show its application to real polarimetric SAR 
data. We also show that using roughness information from the three intensity 
components increases the discriminatory capability of the technique. 

The structure of this paper is as follows. Section |2]presents the new statis- 
tical model for polarimetric SAR data. Section [3] specifies the criterion used 
to determine the transition points and explains the region fitting algorithm, 
including the methodology employed for assessing the precision of the local 
edge detection technique. Sections 3^4] and [4] present the results of applying 
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the technique to real SAR imagery and the conclusions, respectively. Two 
appendices present details about single- channel SAR data (Appendix A) and 
about polarimetric laws, including algorithms for sampling from the random 
variables employed in our description (Appendix B). 



2 Polarimetric SAR Data 

Polarimetric SAR systems use antennae designed to transmit and receive 
electromagnetic waves of a specific polarization, being the two most common 
ones the horizontal linear or H, and vertical linear or V. Due to the possible 
change in polarization of the scattered wave, radar antennae are designed to 
receive the different polarization components simultaneously and, therefore, 
HH, VV, HV and VH data will be available in a full polarimetric system. HV 
and VH channels are identical in an ideal and perfectly calibrated monostatic 
radar system, so one of them will be discarded in the following. 
Then, we define the complex random vector Z as 

Z = [Z HH , Z HV , Zyyf, (1) 

where t V denotes transposition and Z HH , Z HV and Z vv denote the corre- 
sponding components of the backscattered electromagnetic fields, with the 
first subscript indicating the polarization of the transmitted electromagnetic 
field and the second subscript indicating the polarization of the detected 
component of the backscattered electromagnetic field. 

Multi-look signal processing is frequently applied in order to enhance the 
signal-to-noise ratio. We define the multi-look complex matrix of size 

3 x 3 as 

1 n 

ZW = -J^Z(A;)Z* t (A;), (2) 

U k=l 

where V denotes the complex conjugate, n is the number of looks and Z(/c), 
1 < k < n, are random vectors of the form defined in ([I]), so each term of 
summation in ^ is given by 



Z{k)Z*\k) 



Zry* ry ry* ry ry* 

HH^HH ^HH^HV ^HH^VV 

Zry* ry ry* ry ry* 

HV ^ HH ^HV^HV ^HV^VV 

ry ry* ry ry* ry ry* 

^VV^HH % Z W ^VV^VV 



4 



In this Hermitian matrix the diagonal elements are real numbers, while the 
off-diagonal elements are complex numbers with non-null imaginary compo- 
nents. 

We will follow the multiplicative paradigm, so the returned polarimetric 
data will be considered to be the result of the product between backscatter 
variability and polarimetric speckle, given by 



(3) 



where Z = [Zhh, Zhv, Zvv] an d Y = [Yhh, Yhv, Yyvf are complex random 
vectors [29]. The random variable X is scalar, has unitary mean, and models 
backscatter variability due to the heterogeneity of the sensed area. The 
random vector Y represents the polarimetric speckle noise and the mean 
values of each of its components determine the mean values of each of the 
corresponding components of Z. 

The multi-look polarimetric speckle matrix Y^ is defined as 



Zhh 




' Y HH ' 


Zhv 


= Vx 


Yhv 


Zyy 




Yyy 



1 n 

Y (n) = -J2 Y ( k ) Y *\ k )- 



(4) 



k=l 



Then, from g, Q and Q one has that Z< n > = XY {n \ so equation @ can 
be rewritten as 

X n 

Z (n) = _^Y(£;)Y**(A;), 
k=i 



n 



where 



Y(A;)Y* t (fc) 



YhhY hh 
YhvYhh 
YvvYhh 



YhhYhv 
YhvYhv 
YvvYhv 



YhhY vv 
YhvY^v 
YvvYhv 



2.1 A new polarimetric distribution 

If we consider that the components of Y(fc) exhibit a Multivariate Complex 
Gaussian distribution, then riY^ will h ave a Centered Complex Wishart 
distribution, as presented in Appendix B So the density function of Y^ n ^ is 
given by 



/yw (y) 



n \y\ 



exp(-nTr(S Y 1 y)) 



7r 3 r(n)---r(n-2) |S Y | n 



(5) 
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for n > 3 and for Y G C 3x3 , where Tr is the trace, | • | denotes the determinant 
of a matrix and £y is the covariance matrix of Y. 

In order to find the density function of 7S- n \ the following integral has to 
be computed: 

fz(")\x=x( z ) fx(x)dx. (6) 



fzw ( z ) = 
Using (|5| we find that 

fz(")\X=x( z 



X 



'fy(n){x 1 Z). 



(7) 



We consider that the variability of backscatter, modeled by random variable 
X, follows an Inverse Gaussian distribution |39j with unitary mean X ~ 
IG(u, 1) (see Appendix B) with density function 

co (x — 1) 



fx(x) 



2ttx 3 



cxp 



x 



(8) 



where x > and u > is the roughness parameter. This situation, denoted 
X ~ IG(u, 1), is a particular case of the Generalized Inverse Gaussian dis- 
tribution [27] , whose use for backscatter modeling was proposed by Frery et 

al. ng. 

Now, from ^ and (jlj) we obtain 



/z(«) ( z ) 



n 



in— 3 



K, 



3n+l/2( 



u;(2nTr('Sv 1 z) 



! r(n) • • • r(n - 2) |S Y | n ( w (2nTr(E Y 1 z) + 



TV' 



We denote this situation ~ Qp (co, Sy): the Harmonic Polarimetric dis- 
tribution. Appendix A presents the intensity channel version of this and 
other related laws, while Appendix B provides a complete account of the 
polarimetric laws derived from the multiplicative model and the Generalized 
Inverse Gaussian distribution for the backscatter. 

The Bessel function K 3n+ i/ 2 above can be computed using a closed for- 
mula: 

(np + k)\ 



I — np 
y k=0 



(np — k)(2v) k 



with 

v = yj w(2nTr(S Y 1 z) + to), 

alleviating, thus, the numerical issues that the evaluation of this function 
imposes in the general case. 
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2.2 Parameter Estimation 



The estimation of the roughness parameter oj can be done using the first and 
second order moments of the diagonal elements of Z^. The components of 
the main diagonal of TJ^ are given by 

X n 

z\f = -J2 1 1 ^' 2 > with * e ( HH > HV ' vv ^ ' 
fc=i 

where the random variables X and n~ l Y12=i \Yk,i\ 2 are such that X ~ 
IG(u, 1) and n -1 Ylk=i l^k,»| 2 ~ cr |r(n, 2n), where of is the mean of Zj" 
This is equivalent to considering zj^f as the result of the product of a 
IG{uj, of) distributed random variable and a 2n) distributed random 
variable, because a afIG(u, 1) distributed random variable is IG(cu, af) dis- 
tributed. This, implies that Z^f is a Gf{u), erf, n) distributed random variable 
(see Appendix A). Thus, we can estimate these parameters as in the case of 



the univariate intensity data. 

Then, the rth-order moment of the return Z is 




r» ■ 

Now calling rriu = K[Z^} and m 2 j = E[(Z^™' 1 ) 2 ], estimates of u are given by 



1 

Ui -- 



n m,2i ' 

n+1 mfi 



and the estimates of of are given by 

of = m u (9) 

for i G {HH,HV,VV}. 

The value of the parameter u common to the three components will be 
chosen as the one that minimizes the total error e given by 

e = Yl Yl (f z < ( w ' ^ z ki)~ h , 

with i G {HH,HV,VV}, ki G {1, ...,iV}, where iV is the sample size, 
Zi ~ {7p (c^o^n) and /i is a histogram. The parameters that compose the 
correlation matrix are given in table [1} 
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Table 1: Statistical parameters of the correlation matrix. 





HH 


HV 


VV 


HH 




O'HHHV + jbuHHV 


O'HHVV + j^HHW 


HV 






a Hvvv + jbuvvv 


VV 






°vv 




HHVV* -> a HHVV + jb HHVV 
HVVV* -> a H vvv + jbnvv 



HHHV* ->■ a ffray + jb HHHV 



a 



vv 



HV 



HH 



Figure 1: Example of covariance matrix and u parameter estimation for a 
particular region. 



Figure [T] shows the covariance matrix and the images necessary to es- 
timate the statistical parameters of the covariance matrix for a particular 
region, where the arrows relate the images to the corresponding estimated 
parameters. As the covariance matrix is Hermitian, only the upper triangle 
and the diagonal are displayed. 

The off-diagonal elements Z^f in 71^ are given by 

-i n i n 

z if = -E z ^ z h = X -H for * + ^ and *> 1 e i HH > HV > vy } ■ 

k=l k=l 

Due to the independence between X and Y{, with % G {HH, HV, VV}, and 
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taking into account that E [X] = 1, we have that 

= E [ZiZ}} = E [X] E [YiY e *] = E [Y;Y/] . 

As E [YjYg*] = (an + jbit) <7icre, then an and bu can be estimated as 

an+jbu = — , 

where E [ZiZ%] and cjj are estimated from the observed data Z = [Zhh, Zhv, ZwY 
and using equation (|9f , respectively. 

2.3 Parameter interpretation 

One of the most important features of the Q H (both intensity and polari- 
metric) distribution is that the estimated values of the parameter u have 
immediate interpretation in terms of roughness. For values of u near zero, 
the imaged area presents very heterogeneous gray values, as is the case of 
urban areas in polarimetric SAR images. As we move to less heterogeneous 
areas like forests, the value of u grows, reaching its highest values for ho- 
mogeneous areas like pastures and certain types of crops. This is the reason 
why this parameter is regarded to as a roughness or texture measure. 

The parameters of, with i £ {HH, HV,VV}, are the average intensities 
for each polarization. The parameters a^+j b^, with k, £ G {HH, HV, VV}, 
are the correlation coefficients among the three polarimetric images. 

In order to check the capability of the Qp model for describing polarimet- 
ric data, an E-SAR image of WeBling (Bayern, Germany) was used [25] . This 
single look image, shown in Figure |2j was obtained in L band, and it exhibits 
an airport, urban areas, forest and pastures. From this one look complex 
polarimetric image, three intensity images were generated by taking one of 
every three columns (azimuth direction). These images were then averaged 
yielding a three looks intensity image. The parameters were estimated using 
samples from the last three targets, which are shown in Figure [2j 

The estimated covariance matrices for urban, forest and pasture 
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given in (10), (11) and (12), respectively: 



962892 



360932 



19171 - j3579 
56707 

11050 + j3759 
98960 



-154638 + j'19138J 
-5798 + jl6812 
472251 

63896 + j 1581 
6593 + j'6868 
208843 



32556 



556 +J787 
1647 



24046 - j27287 
-146 - j482 
61028 



(10) 

(11) 
(12) 



Table [2] shows the estimated values of u for urban, forest and pasture 
targets, respectively, in HH, HV and VV polarizations. These values are 
comparable in different frequencies when estimated in areas with similar 
textures. Recall from equation (|8j) that to is only allowed positive values; 



c.f. |Appendix B|for more details. 



Table 2: Estimated roughness parameters 



Q 


HH 


HV 


VV 


Average 


Urban 


0.17 


0.94 


0.19 


0.43 


Forest 


10.22 


8.53 


10.55 


9.77 


Pasture 


19.88 


22.54 


18.32 


20.24 



It is noteworthy that the estimated values of the parameters that index 
any polarimetric distribution derived from the multiplicative model do not 
depend, in principle, from the number of looks; the roughness of a target is 
preserved, regardless that measure of the signal-to-noise ratio. 

Figures [3j [4] and [5] show the histograms and estimated densities for differ- 
ent components and zones. It is noticeable that the fit is excellent in urban 
samples and very good in the other cases. 

Figure [6] shows synthetic images with three regions generated with the 
estimated covariance matrices shown above. The background in these figures 
has a 2 HH = a 2 HV = ay V and no correlation among components. 

The estimated values presented in equations (10), (11) and (12), and in 
Table [2] will be also used when assessing the error of the proposed technique. 
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Figure 2: L-band E-SAR image from a zone near the city of Wefiling, Bayern, 
Germany with three regions of interest: Urban (1), Forest (2) and Pasture (3) 

3 Boundary Detection 

In this section we describe an algorithm developed for boundary detection us- 
ing B-Spline deformable contours. Gambini et al. [17] proposed this method 
univariate real SAR imagery, and here we adapt it to polarimetric data. 

B-Splines are a convenient representation of spline functions where the 
curve is specified by a few parameters, the control points; this reduces the 
computational effort to compute it. The order of the polynomial segments 
is chosen arbitrarily, and it relates to the desired smoothness. The B-Spline 
approach allows the local control of the curve by controlling the control 
points individually. The curve lies within the convex hull induced by the 



11 



(a) HH g% (b) HV gg (c) vv g% 

Figure 3: Histograms and estimated densities for urban data 




( a ) hh (b) hv gp (c) vv gp 



Figure 4: Histograms and estimated densities for forest data 

control points. For details of B-Spline representation of contours see the 
works [21 QS1 EE]. 

A more sophisticated approach could be based, for example, on the min- 
imum description length (MDL) principle but the one presented here 
provided excellent results in real SAR applications. 

3.1 Initial Regions 

The procedure begins with a rough segmentation, either manual or automatic 
(with low computational cost), to be refined. 

Let £ be a scene made up by the background B and a region R with 
its boundary dR. We want to find a curve Cb that fits dR in the image. 
We define first an initial search area, which is specified by polygons whose 
vertexes are the control points that generate a B-Spline curve. We developed 
an automatic algorithm for finding initial regions, but it also can be specified 
by the user. 
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(a) HH QU, (b) HV Q§ (c) VV Q§ 



Figure 5: Histograms and estimated densities for pasture data 

If automatic initialization is chosen, the following process is performed. 
The polarimetric image is of the form / : S — > IR 6 , where S — [0, . . . , msu — 
1] x [0, . . . , nsu — 1], m, n, su G N, i.e., it is composed ofmxn blocks of side 
s be . The following data entries have to be specified 

1. A selection criterion Tr C M + which depends on the homogeneity of the 
zone of interest. For example, in order to find urban areas we choose 
T R = [0.1, 1.5); this can be done using natural language and a table 
that converts text into parameter intervals. 

2. A threshold T s that corresponds to the minimum number of blocks 
required for considering a candidate zone as an initial region. This 
specification can also be done in an intuitive and natural manner in 
terms of metric units provided the pixel resolution. 

The parameter oj is estimated for each block S^-, i = 0,...,m, j = 
0, . . . , n, forming an array of size mxnof roughness estimates u(i,j). No- 
tice that each estimate is based on s\ t samples. We opted for su = H, 
after experimenting different values in images with varying complexity; big- 
ger windows provide more precise estimates when acting on areas evenly 
occupied, but will be more prone to mixing samples from different targets. 
If Cj{i,j) G Tr, the block Sij is marked as candidate zone, else it is left 
unmarked. 

If the number of connected blocks of a candidate zone is below Ts, then 
the zone is considered as noise and it is discarded. The initial regions are 
formed by blocks, whose convex hull is then computed. 

In this way, we define an initial search area by means of the automatic 
determination of candidate connected components, which are specified by 
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(a) 



(b) 



(c) 






(d) 



(e) 



(f) 



Figure 6: Synthetic images generated with different covariance matrices: 
a) \HH\ 2 , b) \HV\ 2 , c) \VV\ 2 , d) \HHHV*\, e) \HHVV% f) \HVVV*\ 



polygons whose vertexes are control points that generate a B-Spline curve. 
Alternatively, the user is allowed to manually specify as many as desired 
regions of interest. 

Once the initial search zones are determined, their centroids are calculated 
and the algorithm proceeds with the contour detection. 

3.2 Contour Detection 

If a point belongs to the object boundary, then a sample taken from the 
neighborhood of that point should exhibit a change in the statistical param- 
eters. We consider N segments i — 1, . . . , TV of the form = CPi for 
each candidate area, being C the centroid of the initial region, the extreme 
Pi a point outside of the region and 9 = Z(s^\s^ +1 ^) the angle between 
two consecutive segments, for every i. It is necessary for the centroid C 
to be in the interior of the object whose contour is sought. The points Pj, 
i — 1, . . . , N are arbitrarily chosen with the condition that they are outside 
the object of interest. 
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The segment is an array of m x 6 elements coming from a discretization 
of the straight line on the array of the polarimetric image and is given by: 



,0 



) , 1 < i < A, 



where z£\ k = 1, . . . , m is an array of 6 elements as Figure [7] shows. 




(0 



Figure 7: Data structure for the segment s 



In order to find the transition point on each segment the parameter to 



is estimated as explained in section 2.2 using a rectangle around the segment 
and a sliding window. Then a set of estimates fiW = (u5i, . . . , cD m ) is obtained 
and the biggest variation of ti^ within the array is found, following Blake and 
Isard [2], convolving it with an appropriate mask. The coordinate at which 
cD exhibits the most intense variation is then considered a border point. Once 
the set of border points A = {b%, . . . ,b^} is found, the method builds the 
interpolating B-Spline curve. Algorithm [l] shows a summary of the process 
to find the border points. 



Algorithm 1 Boundary Detection 



1: 



Find or specify initial candidates and use the vertexes of the initial regions 
as control points of the starting boundary. 
Determine a series of radial segments on the image, 
for each segment do 

Generate sets of estimates of the parameter u using the data within a 
rectangular window around the segment. 

Detect the border point by convolving each set with a border detection 
operator. 
6: end for 

7: Return the B-Spline that interpolates the points found. 
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3.3 Error Evaluation 



In this section we present the methodology employed for assessing the preci- 
sion of the edge detection technique previously discussed. 

In Polarimetric SAR images it is very difficult to evaluate the committed 
segmentation error due to the difficulty in establishing the true or correct 
segmentation. Udupa et al. [12] propose a framework for evaluating image 
segmentation algorithms and a systematic way to compare two different al- 
gorithms. In order to obtain accuracy in error evaluation, the authors define 
the surrogate of truth in two possible ways, manual delineation and mathe- 
matical phantoms. The first consists in tracing object boundaries manually, 
while the second consists in creating a set of as realistic as possible simu- 
lated images. In this paper we use the second option, in order to measure 
de local error. We successfully used this methodology for error evaluation in 
univariate amplitude SAR data [T8] . 

A phantom image was used to compute the error in estimating by b the 
true boundary point Pt- This image is a 20 x 100 pixels data set divided 
in halves, and each half is filled with samples from the Qh distribution with 



different parameters using the algorithm presented in B.3 Figure [8] shows 
the one look situation, with u = 10 to the right and u = 1 to the left, both of 
them have the covariance matrix S u , estimated using real data. The figure 
also shows the correct border (the white vertical line at 50) and an estimated 
transition point (the white dot denoted b at 54). The error in this situation 
would be of four pixels. 

Two hundred replications were made for each situation, the transition 
points were estimated and the error was estimated. The distance of these 
points to the true boundary was evaluated and then the array E was defined 
in each situation, as 



E(r) = |50-6(r)|, 1 <r < 200, (13) 

where b(r) is the transition point found in the r-th replication. 

We use relative frequencies in order to estimate the probability of having 
an error smaller than a certain number of pixels. Denote by H(k) the number 
of replications for which the error is smaller than k pixels, then an estimate 
of this probability is given by f(k) = H(k) /200 for k > 0. Algorithm [2] 
illustrates this process. The bigger this probability, the better the algorithm. 
Note that values of k close to are the ones that must be taken into account 
to evaluate the technique under study; gross errors are not interesting. 
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Algorithm 2 Boundary position error estimation 



for each situation do 
for each 1 < r < 200 do 

Simulate a polarimetric SAR image of size 20 x 100, with two regions. 



4: Find the transition point on the major axis of the sample b(r). 

5: Find the distance between the found and the true transition points. 



G 
7 
8 
9 
10 



Update E, as defined in eq. (13) 
end for 
Compute H. 

Return the relative frequencies /. 
end for 



Twelve situations were considered in the study, each one referring to a 
pair of of different type modeled as 

1. E u and oj G {1, 5} besides Ef and oj G {10, 15}; 

2. E u and oj G {1, 5} besides E p and oj G {20, 25}; 

3. Ef and oj G {10, 15} besides E p and oj G {20, 25}. 

The covariance matrices E u , Ef and E p are the ones estimated using real data 



an presented in equations (10), (11) and (12), respectively. 

Figure [9] shows the probability of finding the border point with an error 
lower than the number of pixels indicated on the horizontal axis for each of 
the twelve situations considered. Four curves are shown for each situation: 
three discontinuous, related to single- channel data, and a continuous, which 
presents the results of using the mean of the estimated roughness parameters 
oj on the three channels. 

The first conclusion is that using information from the three channels 
greatly enhances the discriminatory capability of the algorithm, in complete 
accordance with other results in the literature. With the sole exception of sit- 
uation XI, where every technique fails to detect the edge and we observe just 
random fluctuations, continuous lines are above the other ones. Regarding 
the information content of each individual channel, HV polarization outper- 
forms the other two: notice that the curve labeled (c) is in most situations 
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Pi 



50 



UJ 



UJ 



10 



Figure 8: Sample of the phantom image with true and estimated boundary 
points 



above curves (b) and (d). In the following, we will only analyze the results 
obtained with the three channels, i.e., continuous lines labeled (a). 

It is noticeable that urban patches can be easily discriminated from forests 
and from pasture, since all situations but X, XI and XII, where there no urban 
data are used, rapidly rise to values close to 1. The discrimination between 
pasture and forest only is a hard task, c.f. situations X, XI and XII, being 
the first, i.e., £f, u = 10 against S p ,o; = 25 the only feasible among them. 



3.4 Application to real data 



Figure [10] shows the HH band of a 3 looks real E-SAR image showing an 
urban area from the city of Munich, along with a region boundary detected 
using Algorithm [TJ The control parameters specified a medium size homoge- 
neous area and, as can be seen, the technique deals well with both complex 
structures and noisy data. 

Figure [IT] shows the result of applying Algorithm [TJ to the same image, 
but specifying four regions of interest manually. 

Smooth curves were sought in both cases. This was specified setting the 
degree of the polynomials to four. 

Our proposal employs a statistical model with interpretable parameters, 
leading to more information than just the detected edges and the usual ge- 



ometrical features (area, shape etc.). In the case of figure 10 the algorithm 



18 



informs that Q = 17.32; the detected area conforms, then, to our require- 
ment of being homogeneous. The four regions found in figure 11 returned, 
from left to right, the following roughness estimates: Q = 0, 75, u = 18, 94, 
Q = 16, 08 and Q = 12, 88. Those areas are thus labeled as extremely hetero- 
geneous, two homogeneous and one homogeneous tending to heterogeneous, 
respectively. 



4 Conclusions 

Polarimetric SAR imagery segmentation is a very difficult task to solve. In 
this work we described a new approach to region boundary detection in 
polarimetric SAR images using B-Spline deformable contours and local pa- 
rameter estimation. The boundaries of several regions with varying degrees 
of complexity were obtained using our proposal. 

In the first step we either find or specify regions of interest that correspond 
to areas with different degrees of homogeneity, as a coarse first approxima- 
tion. For each region, its boundaries are considered as the initial solution for 
the border detector. Then, the estimated parameter of roughness is calcu- 
lated using two samples: one included in the region and the other out of the 
region and we find the transition point only for the data that are on a set 
surrounding a line segment. All these processes diminish the computational 
cost and improve the performance of the method. 

For each region, the result of the application of this algorithm is a bound- 
ary curve given by an expression in terms of B-Spline functions. The results 
using both simulated and real SAR images are excellent and were obtained 
with an acceptable computational effort. 

In addition, the error in finding edges was defined, and a Monte Carlo 
experiment was used to assess this error in a variety of situations that ap- 
pear in the practice of polarimetric SAR imagery analysis. The information 
content was quantified, in the sense that these results show that using the 
estimator of u computed from the three intensity components consistently 
leads to better results than employing one polarization. 

We did not consider the other nine parameters of the polarimetric model 
for the detection of transition points, given the very good performance ob- 
tained using the estimation of the u alone. 

Future work includes the use of more parameters for finer detail detection, 
improved and robust estimators (as, for instance, the ones presented in PQEl 
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EJ EEU HSl EH S3]) and other types of deformable contour methods based on 
level sets and on stochastic distances [33J . 
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Appendix A Intensity Univariate Distributions 

for SAR Data 

The multi-look return in intensity univariate SAR images can be modeled as 
the product of two independent random variables, one corresponding to the 
backscatter X and other to the speckle noise Y. In this manner Z = X ■ Y 
models the return Z in each pixel under the multiplicative model. For uni- 
variate intensity data, the speckle noise Y is modeled as a T(n, n) distributed 
random variable, where n is the number of looks, while the backscatter 
X is considered to obey a Generalized Inverse Gaussian law, denoted as 
JV- 1 (a,A, 7 ). 

For particular values of the parameters of the N~ x distribution, the 
r(a, A), the r _1 (a,7), and the IG(j, A) (Inverse Gaussian) distributions are 
obtained. These, in turn, give rise to the K, the Q°, and the Q H distributions 
for the return Z, respectively. 

Given the mathematical tractability and descriptive power of the Q° (c.f. 
references (30J I2U [35] ) and the Q H distributions, they represent an attractive 
choice for SAR data modeling. As in this work we will use the Q H we will 
describe it in more detail here. 

The density of the Generalized Inverse Gaussian distribution is given by: 

with parameters 7, A and a in the following parameter space: 

7 > and A > if a < 
7 > and A > if a = 
7 > and A > if a > 0, 
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where 



1a(x) 



1 if x e A 
if x A, 



and where K a denotes the modified Bessel function of the third kind and 
order a, given by 



K a (V ab) 



aW 2 1 



x a x exp 



{ax + bx x ) J . 



The Inverse Gaussian distribution IG(j, A) is obtained when a = —1/2, 



and its density function is given by Eq. (14) 

fx(x) = 



7 
2ttx' 3 



exp \ Yx — r 



R+ {x) , 



(14) 



with A, 7 > 0. The formula for the moments of this distribution is 



E[X r 



7T 



1/2 

7 A ) exp(v / 7^) 



K„ i 



so, the moments of first and second order, and the variance are 
E [X] = ^ E [X 2 ] = 1 + and E [(X - E [X]) 2 ] = 



respectively 

The parameters 7 and A can be used to define a ne w p air of parameters 
to and 77, using w = a/tA and 77 = 
as 



fx(x) 



CUT] 

2ttx 3 



, so equation (|14|) can be rewritten 

2' 



fxp I ^- ) 1 R+ (2) 

2 XT] 



Then, if X ~ JG (tu, 77) it is possible to see that the corresponding moments 
are 

E[X r ] = ^e"r l r Ki(co), 



7T 



so the first and second order moments, and the variance are 

E [X] = 77, , E [X 2 ] = t? 2 ^^, and E [(X - E [X]) 2 ] = t? 2 -, 

" CO ' CO 
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respectively. With this re-parametrization, r/ is the mean value and the 
variance grows as u — > 0. If X ~ IG (u, rj) then X/r] ~ IG (w, 1), and thus 
the density function of the random variable X/rj is given by Eq. (15) 



uj r lo (x — 1 



,2 



( x ) = A/^ exp< ! ~o — i — r 1 ^ ( x ) ( 15 ) 



,r 



In Figure 12 the curves corresponding to this density for rj = 1 and vari- 
ous values of u, are shown. It is noticeable that the variance grows as the 
parameter uj approaches 0. 

Figure 13 exhibits the curves for oj = 1 and various values of 77. Here, the 
curve flattens as the value of 7] grows. 

The density function for the return Z under this model is given by 



with u,r], z > and n > 1. The r-th moment of the Q H distribution is 

Tin) ' 



E gH (Zn=(l) r e«J^K r _ 1/2 (u) 



which is used for parameter estimation. The modified Bessel function of the 
third kind and order v, whose integral representation is, according to [24], 
given by: 

/■oo 

K u {z) = I exp{— z cosh(t)} cosh(z/t)c?t 
Jo 

is here denoted K v . Numerical problems arise when computing this function 
(c.f. [22]), but the distribution used in this paper circumvents this issue as 
will be seen in the next section. 



Appendix B Polarimetric Laws under the Mul- 
tiplicative Model 

This section presents the distributions for polarimetric SAR data: the Com- 
plex Multivariate Gaussian distribution and the Centered Wishart distribu- 
tion (see [2Ql [2H [HI E]). This last one is the most frequently used model 
for the return coming from homogeneous areas, and serves as a basis for the 
return from heterogeneous and very heterogeneous areas [34] . 
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B.l Complex Multivariate Gaussian distribution 



Let T = [ti, . . . ,t m ] be a random m-dimensional complex vector with Gaus- 
sian Complex Multivariate i.i.d components. Each component has the form 
tk = Rk + Hk where Rk and Ik are the real and imaginary parts, respectively, 
and they are real random variables. Then, we can define a 2m-dimensional 
Gaussian Multivariate vector H = [R±, I±, . . . , R m , I m ]. Its covariance ma- 
trix Eh is a block matrix with blocks given by (En)fc/ of size 2x2 with 
k, £ — 1, . . . , m, as follows 



(Rk 



Hn k )(Re - (J>r. 
1 



(Rk 
(h 



m k )(h 
Hi k )(h 



Hit) 



2 





bke 



1 

—b k e 
&ke 



if k = £ 



if k^£ 



where o- k /V2 are the standard deviation of each component of the random 
vector H, and bki are the correlation coefficients. 

If the random vector H is multivariate normal distributed, then the cor- 
responding random vector T follows a Complex Multivariate Gaussian dis- 
tribution, which denoted by T ~ He (a*t, S t ). The density function is given 
by 



h(T) 



1 



exp (- (T - /! T )** S t x (T - /i T )) , 



where /i T is the mean value and the covariance matrix 

if k = £ 
if k^£, 



T)k > e \ (a k e + jbki) (JkOe 



(16) 



with k,£ — 1, . . . , m. 

The Multivariate Complex Gaussian distribution is the base of the Cen- 
tered Complex Wishart distribution which is the return model corresponding 
to homogeneous areas. 



B.2 Centered Complex Wishart Distribution 

The Centered Complex Wishart Distribution is used to model the speckle 
noise of polarimetric data, the n looks are considered as n random vectors 
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T (1) , . . . , T (n), i.i.d, with T (k) ~ J\fc (0, St) whose dimension is m, m < 

71. 

The random matrix W of m x m is defined as: 

n 

w = J2 T ( k ) T ( k T t - ( 17 ) 

fc=i 

Then, the joint distribution of the m x m elements of the W matrix 
is the Centered Complex Wishart distribution [19] and is denoted as W ~ 
W (ET,n), and the parameter n indicates the degrees of freedom. 

The density function of the random matrix W is given by: 



/w (W) = z^ivaw-Y 1 ^ ~ , TYTvCT" ex P H r £t ^)) 





W 


n— m 


7r m(m-l)/2 r ( n ) . 


■ ■ r (n — m + 1) St 


n 



for n > m and for all G C mxm . 

In polarimetric SAR images analysis this distribution is used to describe 
homogeneous areas. The parameters, that characterize each different region 



on the image, are the matrix values given by the equation (16). 



In order to estimate the distribution parameters, we use the principal 



diagonal components of the W, defined in equation (17) and given by: 



W hl = J2\ T *( k )\ 2 > ie{l,...,m} 



i=i 



It holds that W ¥ ~ nafT{n, 2n), with af = E [\Ti{k)\ 2 ] for all k = 1, . . . , n. 

If a random variable W has a T(n, 2n) distribution, then E(W) = 1, and 
^(W^i) = naf for every i. 

From the moments method the estimator of a,- results 



n 



If the areas have different degrees of heterogeneity, it is necessary to gener- 
alize this model introducing the possibility of having variable characteristics 
instead constant features. To this end, Yueh et al. [45J proposed the K polari- 
metric distribution. In the practice, we find extremely heterogeneous data, 
for this reason a more flexible and tractable distribution was proposed, called 
the Harmonic polarimetric distribution and denoted as Qp . 
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B.3 Drawing outcomes from the Q distribution 

The Q H distribution belongs to the multiplicative model, and it describes 
the law that governs the product Z = XY, where the independent ran- 
dom variables X and Y follow the Inverse Gaussian and Complex Wishart 
distributions, respectively. 

Outcomes from the Complex Wishart distribution can be obtained from 
transformations of Complex Gaussian distributed outcomes, while outcomes 
from the Inverse Gaussian distribution are obtained by an acceptance- rejection 
technique. These procedures are described in the following. 

B.3.1 Multivariate Gaussian Random variable generation 

If a random vector Y = [Yj., . . . , Y p f is Multivariate Normal distributed then, 
its density function is given by 

/y (y) = (27r) n/ "\j: Y \ 1/2 exp (~ 2 1 (y " ^ Sy ' (y " ^ y) ) ' 

where y = [yi, . . . , y p ]*, /Iy is the mean vector of Y, and S Y is the covariance 
matrix, which is a symmetric positive definite matrix. 

The elements of £y are given by s^- = p^a^aj, 1 < i,j < p, where <7j 
is the standard deviation of the random variable Yi, and p^ is the correlation 
coefficient between Yi and Yj. It is possible to verify that \pij\ < 1, that 
Pij = pji and that pa — 1 for all 1 < i < p. 

Let $y be the p x p-dimension matrix whose columns are the normalized 
eigenvectors of the matrix Sy and let Ay be the diagonal matrix with the p 
eigenvalues of the Sy in the diagonal elements, then we have Sy'I'y = $yAy- 

In order to generate multivariate normal random values Y with mean 
value py and covariance matrix Sy, a set of decorrelated zero- mean normal 
values W (£w = I and pw = 0), are generated. They are the transformed 
by Y = $ Y Ay /2 W + p Y . 

B.3. 2 Inverse Gaussian distribution generation 

Algorithm [3] shows how to generate samples from the IG(ui,r]) distribu- 
tion jZj. 
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Algorithm 3 Inverse Gaussian distribution generation. 
1: Generate t, sample of the random variable T ~ A/"(0, 1) 
2: Calculate v = rj + £ - ^y/i (Au + t) 
3: Generate u, sample of the random variable U ~ W(o,i) 
4: if u > 1/(1 + U77) then 
5: return (r^w) 1 
6: else 
7: return v 
8: end if 
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Probabilities 




Figure 10: Real polarimetric three looks E-SAR image and boundary detec- 
tion with automatic initialization 
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Figure 11: Result of applying Algorithm [T] to a 3-looks polarimetric image 
with four initial regions manually specified 
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1.4: 




Figure 12: Density of the IG(x,oj,rj) distribution, for 77 = 1 and oj = 1 (solid), 
oj = y/2 (dashed), oj = 2 (dotted), oj = 3 (dot-dash), w = 4 (dot-dot-dash), w = 5 
(solid) and oj = 6 (dashed). 
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Figure 13: Density of the IG (x,uo,rj) distribution, for cj = 1 and rj = 0.5 (solid), 
r] = 1 (dashed), i] = 2 (dotted) y rj = 3 (dot-dash). 
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